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ABSTRACT 

We present the results of N-body simulations of dissipationless galaxy merging in 
Modified Newtonian Dynamics (MOND). For comparison, we also studied Newtonian 
merging between galaxies embedded in dark matter halos, with internal dynamics 
equivalent to the MOND systems. We found that the merging timescales are signifi- 
cantly longer in MOND than in Newtonian gravity with dark matter, suggesting that 
observational evidence of rapid merging could be difficult to explain in MOND. How- 
ever, when two galaxies eventually merge, the MOND merging end-product is hardly 
distinguishable from the final stellar distribution of an equivalent Newtonian merger 
with dark matter. 

Key words: gravitation — stellar dynamics — galaxies: kinematics and dynamics 



1 INTRODUCTION 

Given the remarkable ability of MOND to reproduce the 
kinematics of galaxies (e.g., Milgrom 2002; Sanders & Mc- 
Gaugh 2002) and its increased interest due to the pos- 
sibility of a relativistic formulation (Bekenstein 2004), it 
is natural to look for tests able to discriminate between 
MOND and Newtonian gravity with dark matter (DM) . De- 
spite numerous attempts, no clear-cut cases have been found 
so far (see Bekenstein 2006), the main reason being that 
MOND is a non-linear theory, and this makes the study 
of systems deviating significantly from spherical symmetry 
(Brada & Milgrom 1995; Ciotti, Londrillo & Nipoti 2006) 
and the performance of N-body simulations more difficult 
in MOND than in Newtonian dynamics. In particular, very 
few N-body simulations of MOND systems have been per- 
formed so far, exploring stability of disk galaxies (Brada & 
Milgrom 1999; Tiret & Combes 2007), external field effect 
(Brada & Milgrom 2000), dissipationless collapse and phase 
mixing (Nipoti, Londrillo & Ciotti 2007a, hereafter NLC07; 
Ciotti, Nipoti & Londrillo 2007). In addition, the study of 
structure formation in MOND using N-body cosmological 
simulations is still limited to a couple of preliminary explo- 
rations (Nusser 2002; Knebe & Gibson 2004). 

Observations leave no doubt that galaxy merging occurs 
(e.g. Arp 1966; Schweizer 1982), and it is also known that 
Newtonian gravity can account in detail for such a process 
(Toomre & Toomre 1972). It is then natural to study galaxy 
merging in MOND. In fact, there are reasons to expect that 
galaxy merging is less effective in MOND than in Newtonian 
gravity: in MOND galaxies are expected to collide at high 
speed, and there are no DM halos to absorb orbital energy 
and angular momentum (Binney 2004; Sellwood 2004); in 



addition, it has been recently shown that violent relaxation 
and phase mixing are slower in MOND (NLC07; Ciotti et 
al. 2007). 

Taking advantage of our recently developed MOND N- 
body code, in this Letter we present the results of N-body 
simulations of galaxy merging in MOND, focusing for sim- 
plicity on the case of dissipationless merging between equal- 
mass spherical galaxies. For comparison, we also consider the 
merging of structurally identical purely baryonic Newtonian 
systems, and the merging of equivalent Newtonian systems 
(i.e., Newtonian models with the same baryonic distribution 
as the MOND systems, embedded in DM halos such that 
their internal dynamics matches the corresponding MOND 
cases; see Milgrom 2001; Nipoti et al. 2007b). 



2 THE NUMERICAL SIMULATIONS 

We consider MOND in Bekenstein & Milgrom's (1984) for- 



mulation, in which the Poisson equation V - 
substituted by the non-relativistic field equation 



IV0II 
ao 



4nGp. 



AnGp is 



(1) 



In the equation above ||...|| is the standard Euclidean norm, 
4> and (j) are, respectively, the MOND and Newtonian grav- 
itational potentials produced by p, and for finite mass sys- 
tems V4> — > for ||x|j oo. The function /i(j/) is not con- 
strained by the theory except that it must run smoothly 
from p,{y) ^ y a.t y <^ 1 (in the so-called deep-MOND 
regime) to p,{y) ~ 1 at j/ ^ 1, with a dividing acceleration 
scale ao — 1.2 x 1 0~^°m s~^, and in the present work we 
adopt p(y) = y/^Jl + (Milgrom 1983). From the Pois- 



2 Nipoti et al. 



Table 1. Parameters of the simulations and properties of the end-products. 



Name 


Gravity 




K 




bo /do 








c/a 


b/a 


rm/rM.o 


o"v/o"v,o 


7 


(m) 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 




(7) 


(8) 


(9) 


(10) 


(11) 


(12) 


(13) 


(14) 


Mlh 


MOND 





1 


0.958 





2 


X lO" 





0.52 


0.54 


1.6 


1.09 


1.4 ± 0.3 


4.4 ±0.3 


Elh 


Newton 


30 


1 


1.017 





2 


X lO'"^ 


2 X 10^ 


0.48 


0.50 


1.2 


1.22 


1.0 ±0.2 


3.4 ±0.3 


Mlo 


MOND 





1 


0.958 


0.5 


2 


X 10^* 





0.51 


0.61 


1.9 


1.09 


2.0 ±0.3 


6.3 ±0.3 


Elo 


Newton 


30 


1 


1.017 


0.5 


2 


X 10^ 


2 X 10^ 


0.54 


0.60 


1.2 


1.23 


1.3 ±0.2 


4.1 ±0.2 


M25h 


MOND 





25 


0.428 





2 


X 10"^ 





0.54 


0.73 


2.1 


1.16 


1.8 ±0.3 


4.2 ±0.3 


E25h 


Newton 


5 


25 


0.447 





2 


X 10^ 


1 X 10« 


0.52 


0.57 


1.5 


1.27 


1.3 ±0.3 


3.5 ±0.3 


M25o 


MOND 





25 


0.428 


0.5 


2 


X 10'^ 





0.59 


0.79 


2.2 


1.16 


2.0 ±0.3 


6.3 ±0.5 


E25o 


Newton 


5 


25 


0.447 


0.5 


2 


X 10^ 


1 X 10*^ 


0.62 


0.83 


1.3 


1.29 


1.7 ± 0.3 


4.5 ±0.3 


NOh 


Newton 







0.183 





2 


X 10** 





0.68 


0.71 


2.0 


1.07 


2.0 ±0.1 


5.0 ±0.3 


NOo 


Newton 







0.183 


0.5 


2 


X lO*^ 





0.64 


0.89 


1.6 


1.07 


2.1 ± 0.2 


5.5 ±0.3 



(1); name of the simulation. (2): gravity law. (3): DM to baryonic mass ratio. (4): internal acceleration ratio. (5): normalised relative 
speed of the barycentres at i = 0. (6): normalised impact parameter at t = (in all cases do = 40rt). (7); total number of stellar 
particles. (8); total number of DM particles. (9): end-product minor-to-major axis ratio. (10); end-product intermediate-to-major axis 
ratios. (11): final-to-initial half-mass radius ratio. (12): final-to-initial virial velocity dispersion ratio. (13): best-fit inner logarithmic 
slope of the final angle-averaged density profile. (14): best-fit Sersic index of the final projected density profile. 



son equation and equation ([1} it follows that the MOND 
(g = —Vcj)) and Newtonian (g^ = — V<^^) gravitational 
fields are related by jj.{g/ao)g — + S, where g = ||g||, 
and S is a solenoidal field dependent on the specific p con- 
sidered: in general one cannot impose S = 0, thus the use of 
standard Poisson solvers to develop MOND N-body codes is 
not possible, and equation ^ must be solved at each time 
step (Brada & Milgrom 1999; NLC07). 



2.1 Initial conditions and the code 

The baryonic component of the initial conditions of all the 
simulations presented in this paper consists of two identical 
galaxy models with stellar density distribution 



p*(r) = 



A/* 



2tv r{r + r*Y 



(2) 



where M, is the total stellar mass and is the core ra- 
dius (Hernquist 1990). To each MOND model with po- 
tential (j) corresponds an equivalent Newtonian model with 
(j)^ = 0, thus having a DM halo with density pdm(t') = 
V'^0(r)/47rG — Pt{r). In principle, such a DM halo would 
have infinite mass, so we truncate it at r ~ 30r,. For com- 
pleteness, we also ran simulations of Newtonian merging be- 
tween purely baryonic systems with the same stellar density 
distribution and no DM halo. 

The particles of the stellar component are distributed 
with the standard rejection technique applied to the phase- 
space distribution function (DF), restricting for simplicity 
to the fully isotropic case. In the purely baryonic Newtonian 
case the DF is known explicitly (Hernquist 1990), while in 
MOND the corresponding DF is obtained numerically with 
an Eddington inversion (e.g. Binney & Tremaine 1987) 



fM{E) = 



d 



VSn^ dE 



dp, d<j) 



(3) 



where the upper integration limit reflects the far-field loga- 
rithmic behaviour of the MOND potential (see also Angus, 
Fameiey & Zhao 2006). Finally, in the equivalent Newto- 
nian models the stellar particles are distributed by using 



their numerical two-component isotropic DF. However, this 
is not possible in general for the equivalent DM halo parti- 
cles, because for systems with sufficiently high stellar surface 
density the resulting halo presents a central hole, and so it 
cannot be derived from an everywhere positive, isotropic 
DFQ (Ciotti & Pellegrini 1992). Thus, the initial DM parti- 
cle velocities are extracted from a Maxwellian distribution 
with local velocity dispersion satisfying the isotropic two- 
component Jeans' equations. We verified that the resulting 
models are in approximate equilibrium by evolving them in 
isolation for several dynamical times. 

We consider both head-on and off-centre encounters. 
In the head-on cases (impact parameter bo = 0) the two 
galaxies are released at t = with barycentric distance do = 
40r, , and with the relative speed vo that they would have if 
they started at rest at drcst = 60r, . Thus, in the Newtonian 
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4G(M. + Md 
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d^cst 



while in the MOND cases 



■"0 



^-0.8V8GM,ao ln%i. 



do 



(4) 



(5) 



where we have used the approximate expression of the force 
between two particles in deep-MOND regime (Milgrom 1986; 
Milgrom 1994). In the off-centre cases, do and t^o are the 
same as in the corresponding head-on cases, but the relative 
velocity is oriented so that the impact parameter 6o — do/2. 

The physical scales of the problem are introduced as fol- 
lows. First of all, we identify each MOND initial condition 
by fixing a value for the dimensionless internal acceleration 
parameter k = GMt/{aorl), so M» and r« are not indepen- 
dent quantities: in physical units, r, ~ 3.4k~^/^ Af^^p kpc, 
where Af«,io = Afj/lO^^Af©. The time and velocity units are 
t, = ^rl/GKL ~ 29.7k~^/*A'4\/io Myr, and = r-^/U ~ 



^ This result shows that it is important to check the positivity 
of the DF (and not only that of pdm)i when studying Newtonian 
systems with DM equivalent to MOND models. 
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Figure 1. Time evolution of the virial ratio (top panels) and of 
the barycentric relative speed (bottom panels) for the two sets of 
simulations with k = 1; head-on (left panels) and off-centre (right 
panels). Blue, red, and green curves refer to MOND, equivalent 
Newtonian, and purely baryonic Newtonian simulations, respec- 
tively. For the scaling of time and velocity, see Section 12.11 



112«:i/''M,^''iokms-^ (see NLC07 for a more detailed dis- 
cussion of the normalisations) . The simulations are evolved 
up to t = 400t* [k = 1 cases) or t = 500t* (k = 25 and 
purely baryonic Newtonian cases), which amount to several 
gigayears in physical units for galaxy masses in the observed 
range. 

Our MOND N-body code (NLC07) is a parallel three- 
dimensional particle-mesh code that can be used to run 
MOND as well as Newtonian simulations. The code is based 
on a grid in spherical coordinates, on which the MOND po- 
tential is computed by solving exactly the field equation ([T]) 
with the iterative potential solver (based on spectral meth- 
ods) described in Ciotti et al. (2006). Particle- mesh inter- 
polations are obtained with a quadratic spline in each coor- 
dinate, while time stepping is given by a classical leap-frog 
scheme. The time step is the same for all particles and is 
allowed to vary adaptively in time. Given the spherical ge- 
ometry of the grid, the code is not is not best-suited to run 
merging simulations: in order to obviate this difficulty, we 
used a time-adaptive grid with a much larger number of grid 
points (128'^) than in the collapse simulations of NLC07. 
With this resolution, we obtained excellent agreement be- 
tween Newtonian merging simulations run with the MOND 
code and simulations (starting from the same initial condi- 
tions) run with our FVFPS treecode (Fortran Version of a 
Fast Poisson Solver; Londrillo, Nipoti & Ciotti 2003; Nipoti, 
Londrillo & Ciotti 2003) . In summary, all the presented sim- 
ulations were run with the MOND code, and the Newtonian 
simulations also with the FVFPS code. The properties of 
the simulations, including the total number of stellar (Af, ) 
and DM (A'^dm) particles used, are summarised in Table 1. 



Figure 2. The same as Fig. [T] but for the two sets of simulations 
with K = 25: head-on (left panels) and off-centre (right panels). 

3 RESULTS 

3.1 Merging dynamics and timescales 

We present now the results of four sets of merging simula- 
tions, each characterised by a combination of the value of 
the internal acceleration parameter (k = 1 or k = 25) and 
of the impact parameter (bo = or bo = rfo/2). Each set 
comprises three simulations: MOND, equivalent Newtonian 
(A/dm = 30M, for re = 1, and Mdm = 5M. for k = 25) 
and purely baryonic Newtonian. In the top panels of Figs.[T] 
and [2] we show the time evolution of the virial ratio 2K/\W\ 
(where K is the total kinetic energy and W is the trace 
of the Chandrasekhar potential energy tensor) and, in the 
bottom panels, the time evolution of the relative speed Vrd 
of the barycentres of the two galaxies: note that the time 
and velocity units are the same for all simulations. Peaks in 
27^/114^1 and in v^^i correspond to close encounters between 
the two systems, while minima of 2K/\W\ and v-^d occur 
when the separation is maximum. At the end of all the pre- 
sented simulations 2K/\W\ ~ 1 and Wroi ~ 0, indicating that 
the two systems merged, forming a single virialised object. 

Let us focus first on the case k = 1, in which the ini- 
tial galaxies, having internal accelerations everywhere lower 
than ao, are in deep-MOND regime. As can be seen from 
Fig. [TJ in both the head-on and the off-centre cases, the 
merging timescale is apparently longer in MOND (blue 
curves) than in the equivalent Newtonian simulations (red 
curves) . In MOND the two galaxies experience several close 
encounters before merging, while in the equivalent Newto- 
nian cases they merge quickly after the first close passage. 
The behaviour of both MOND and purely baryonic New- 
tonian cases (green curves) is very sensitive to whether the 
orbit is head-on or off-centre: the merging timescale in sim- 
ulation NOo is almost a factor of two longer than in simula- 
tion NOh, and also simulation Mlo takes significantly longer 
to virialise than simulation Mlh. In contrast, due to the 
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presence of DM halos, the merging timescale is as short in 
the ofT-centre as in head-on the equivalent Newtonian cases. 
As expected, the relative speed during the first close en- 
counter is significantly higher in MOND simulations than in 
the purely baryonic ones. On the other hand, the equivalent 
Newtonian models collide at higher speed than their MOND 
counterparts. We note that this last result depends on the 
specific choice of Mom and drost appearing in equations (|4]) 
and ([S]): provided that do is large enough MOND would have 
no problem in attaining arbitrarily high collision speeds (see 
also Angus & McGaugh 2007 for a discussion of the collision 
speed of galaxy clusters in MOND). 

The case k = 25 (Fig. [2}, in which the initial galaxy 
models have internal accelerations >ao for r<5rt, confirms 
the same trend as the a — 1 case, with merging taking longer 
in MOND than in equivalent Newtonian models (by a fac- 
tor of ~ 2 in the head-on case, and by a factor of ~ 4 in 
the off-centre case). The ac = 25 simulations are interesting 
also because they show how the merging process in MOND 
is very difi'erent from that of purely baryonic Newtonian 
merging, even when the MOND galaxies are internally in 
Newtonian regime. This can be easily seen in the head-on 
simulation M25h, in which the dynamics is almost Newto- 
nian when the two galaxies interpenetrate, but the collision 
speed is higher than in the purely baryonic Newtonian case, 
being determined by the long-range deep-MOND interaction 
between the two galaxies. 

We note that the value of k contains information 
only on how the initial internal accelerations compare with 
oq. a model initially characterised by accelerations every- 
where weaker than ao can produce accelerations significantly 
stronger than ao during its dynamical evolution. This be- 
haviour was observed by NLC07 in MOND dissipationless 
collapse simulations (see also Nusser & Pointecouteau 2006, 
who studied spherically symmetric MOND gaseous collapses 
in a cosmological context). However, this is not necessarily 
the case in MOND merging, because during the dissipation- 
less merging process the density does not increase as much 
as in a collapse. To quantify this effect we computed at each 
time step the fraction of particles with acceleration stronger 
than ao. In simulations M25h and M25o this fraction is ini- 
tially ~ 0.3, has a peak up to ~ 0.5 during the first close pas- 
sage, and is again ~ 0.3 in the end-products. On the other 
hand, it turns out that in simulations Mlh and Mlo this 
fraction is <0.02 throughout the entire simulation; in other 
words our n — 1 simulations are in deep-MOND regime at 
all times. 

From an observational point of view, simulations with 
K = 1 can represent merging between two dwarf spheroidal 
galaxies with M, = 10^ (and effective radius Rc ~ 
0.2 kpc, so the merging timescale would be <1.8Gyr), while 
simulations with k — 25 can represent merging between 
two luminous elliptical galaxies with ~ 10" Mq (and 
effective radius -Ro ~ 4 kpc, so the merging timescale would 
be <2.1Gyr). Thus, restricting to the presented cases, one 
could be tempted to conclude that galaxies in MOND can 
merge in a timescale significantly shorter than the Hubble 
time. However, we stress that such a conclusion is wrong 
in a general sense, because we reported only cases with or- 
bital energies corresponding to two galaxies at rest when at 
relatively small distance (drost — 25rM,o, where tm.o is the 
half-mass radius of the initial stellar distributions). Given 



the logarithmic nature of the MOND gravitational poten- 
tial (see equation [5]) , choosing larger values of drcst has the 
effect of boosting the encounter relative speed (making the 
merging process difficult), while it barely affects the en- 
counter speed in Newtonian gravity. In fact, we explored 
several other cases of MOND encounters, with larger drost 
(and then higher vq), but we had to stop the simulations, 
because the two galaxies after the first close passage reach 
relative distances significantly larger than do , making the re- 
quired computational time exceedingly long, revealing virial- 
isation times even longer than the age of the Universe (note 
that the MOND simulation in the right panels of Fig. [T] is 
already dangerously long). Summarising, we presented here 
only simulations of encounters relatively favourable to merg- 
ing in MOND, and yet these mergings were found to be less 
effective than in Newtonian gravity with DM. 

3.2 Merging end-products 

We define merging end-products the systems comprising the 
bound stellar particles at the end of the simulation. In the 
Newtonian simulations we found that <4 per cent of the 
stellar particles escaped, while there cannot be escapers in 
the MOND cases. Using the same procedure as in NLC07, 
we determined the following properties of the end-products: 
the axis ratios c/a and b/a of the inertia ellipsoid, the angle- 
averaged half-mass radius tm, the virial velocity dispersion 
(TV, the inner slope 7 of the 7-model (Dehnen 1993; Tremaine 
et al. 1994) that best fits the final angle-averaged density 
profile (over the radial range 0.1 ^r/rM ^ 10), and, for the 
three principal axis projections, the circularised effective ra- 
dius Re and the index m of the Sersic (1968) law that best 
fits the circularised projected density profile over the radial 
range 0.1 < R/Rc 10 (see Table 1, where (m) is the average 
of the values of m obtained for the three projections). 

The structural and kinematic properties of the MOND 
end-products are not significantly different from those of 
their Newtonian equivalent counterparts: for instance, the fi- 
nal axis ratios are roughly the same in corresponding MOND 
and equivalent Newtonian simulations (see Table 1). The 
end-products of simulations Elh and Elo are DM dominated 
at all radii, and similarly the end-products of the k = 1 
MOND mergers are everywhere in deep-MOND regime (so 
they would appear as DM dominated at all radii if inter- 
preted in Newtonian gravity). On the other hand, the n — 25 
MOND end-products would appear in Newtonian gravity as 
baryon-dominated in the inner regions (r/rM^l.2— 1.3) and 
DM dominated at larger radii; the corresponding equiva- 
lent Newtonian end-products are baryon dominated at radii 
r/rM^O.4 — 0.5. In general both MOND and equivalent end- 
products have rather flat intrinsic and projected velocity 
dispersion profiles. The MOND final density profiles tend to 
be steeper (7 = 1.4 — 2.0, (m) — 4.2 — 6.3) than those of the 
equivalent Newtonian cases (7 = 1.0 — 1.7, (m) = 3.4 — 4.5), 
but there is not a dichotomy between the two families. 

An interesting point (in the context of the galaxy scal- 
ing relations) is how the final virial velocity dispersion ctv 
and half-mass radius tm compare with the corresponding 
quantities in the initial systems ctv.o and rM,o (Nipoti et 
al. 2003). MOND mergers have larger vm and lower av than 
the corresponding equivalent Newtonian mergers. We also 
note that in the Newtonian cases here presented the ratio 
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o"v/o"v,o tend to be larger (and tm/j'm.o smaller) than in 
similar cases explored in Nipoti et al. (2003): this is ex- 
pected, because here we consider elliptic orbits while Nipoti 
et al. (2003) considered parabolic orbits. 



4 DISCUSSION AND CONCLUSIONS 

The main result of the present work is that galaxy merg- 
ing is much less effective in MOND than in Newtonian dy- 
namics with DM. In addition, the derived MOND merg- 
ing timescales must be considered only lower limits, because 
rather specific orbital properties are required in MOND in 
order to have galaxy mergers on timescales shorter than the 
age of the Universe. In general, repeated high speed galaxy 
encounters should be a common feature of galaxy interac- 
tions in MOND, while any observational evidence of rapid 
merging after the first close passage should be regarded as an 
indication of the presence of DM halos. Remarkably, when 
the orbital parameters are favourable and two galaxies even- 
tually merge in MOND, the merging end-product is hardly 
distinguishable from the final stellar distribution of an equiv- 
alent Newtonian merger with DM. 

Thus, the very observation of galaxy mergers appears 
to favour the DM scenario with respect to the MOND 
hypothesis. Additional constraints for galaxy merging in 
MOND could be also given by specific dynamical features 
in galaxy interactions that have extensively studied and ex- 
plained in the context of Newtonian gravity (e.g. Binney 
& Tremaine 1987), such as the tidal tails observed around 
interacting disk galaxies as the "Antennae" pair of galax- 
ies NGC4038/NGC4039 (Toomre & Toomre 1972), and the 
surface brightness ripples observed in the outskirts of lumi- 
nous elliptical galaxies as NGC3923 (Quinn 1984). 

The result that merging is less effective in MOND than 
in a DM scenario appears consistent with our previous find- 
ings that phase-mixing and violent relaxation are slower 
in MOND than in Newtonian gravity (NLC07; Ciotti et 
al. 2007). The merging process is intimately related also to 
dynamical friction, so our simulations might be interpreted 
as an indication that dynamical friction is less effective in 
MOND than in Newtonian gravity with DM, in contrast 
with the analytical estimates of Ciotti & Binney (2004) for 
the case of a particle moving in a homogeneous medium. 
However, the complexity of the merging process prevents us 
from drawing firm conclusions on this issue, and we plan 
to realise ad hoc numerical experiments to explore in detail 
dynamical friction in MOND. 

We must also recall that we explored only very sim- 
ple cases of galaxy merging in MOND: in particular, we 
only considered equal-mass dissipationless merging between 
spherical systems, while dissipative processes in the merging 
of gas-rich galaxies might be effective in making the merg- 
ing timescales shorter. Another possible caveat is that, given 
the long-range nature of MOND gravity, the restriction to an 
isolated pair of galaxies might not be as justified as in New- 
tonian gravity, and the next step to address this point would 
be the study of galaxy merging in MOND in a cosmologi- 
cal context. A valuable contribution in this direction would 
be the performance of cosmological simulations of structure 
formation, based on a self-consistent relativistic formulation 
of MOND such as Bekenstein's (2004) TeVeS. 
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